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Abstract 

We consider the use of a multigrid method with central differencing to solve the 
Navier-Stokes equations for high-speed flows. The time-dependent form of the equations 
is integrated with a Runge-Kutta scheme accelerated by local time stepping and variable 
coefficient implicit residual smoothing. Of particular importance are the details of the 
numerical dissipation formulation, especially the switch between the second and fourth 
difference terms. Solutions are given for two-dimensional laminar flow over a circular 
cylinder and a 15 degree compression ramp. 
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1980’s a wide variety of numerical schemes were investigated for solving 
the eX and NaL-Stokes equations. Multistage time-stepping schemes with centra 
differencing and multigrid acceleration [1, 2, 3] were demonstrated to be quite effect 
Xcomoutfng subsonic and transonic flows over aerodynamic components and config- 
urations. With the recent resurgence of interest in high-speed flig t ve ic es, we now 
need to construct versatile algorithms for hypersonic flow. One must I keep in mind t hat 
hypersonic flows represent a formidable challenge for any flow solver In paiticuiar 
strong tick and expansion waves can occur in the flow field, and hey can interact 

with each other and with shear layers (i.e., boundary layers, jets, ^“numericaTinte 
nonlinear behavior and interactions can easily cause divergence o y 
rtion procedure. This is especially true during the initial phase of a calculation with a 
time dependent method. So even the most successful algorithms of the last decade may 
require significant modifications to be effective for hypersonic flows. 

An inftial effort [4] to apply a central-difference multignd algorithm to high-speed 
flows rX^ in numerical dfltalUe. that prevented the calculation „ two- ,_a 
flows (i e , blunt body and wedge type) with a Mach number higher than about 7. I 
order to compute such flows a low Courant-Friedrichs-Lewy (CFL) number was I - 
Thus four and five stage schemes were not practical, since there ,s substantial dote , 0 - 
ration in the high frequency damping of the scheme due to the arge reduction m the 
CFL number The CFL restriction reduced the potential of the scheme as a viscous flow 
solver. More recently an algorithm utilising a semicoarsening technique, a symmetlic 
TVD formulation, and a three stage Runge-Kutta scheme [5] was propose an 
to compute high Reynolds number (laminar) Mach 10 flow over an airfo. at 10 degrees 
angle Jf attack. A good resolution of the bow shock wave and a reasonable convergent 
rate were obtained. The method of semicoarsening considered required muc 
complicated cycle strategy than that employed with standard multigrid me ° ■ 

addition, it appears to be somewhat cumbersome to implement in three dimensions 

It is our contention that standard multigrid techniques can be used in conjunction 
with central differencing to compute hypersonic flows effectively. To achieve sue success 
with these techniques one needs to give appropriate attention to both the advecti 
and the dissipative processes of hyperbolic multigrid. The advect, on process provides 
a mechanism by which long wave disturbances can be rapidly expelled^ In a multignd 
method with a time-dependent iterative procedure, efficiency is m part derived fiom 
larger time steps allowed on coarser meshes. Hence, it is important that the driving 
scheme of the multigrid method use large time steps. The dissipative process is essen ia 
in smoothing short wave disturbances. With this process the multignd efficiency is based 

on principles similar to that for elliptic equations. 

For hypersonic flows one encounters an additional consideration regarding the damp- 
ing of the short waves. As the Mach number increases the jumps across shocks become 
larger and it becomes more difficult to eliminate these high frequency oscillations. I hus 
a considerable part of the following discussion will concentrate on the smoothing al- 
gorithm. The fundamental features of the multigrid process (i.e.. Full Approximation 
Storage scheme, grid transfer operators, fixed cycle strategy) are fairly standard Other 
aspects, such as type of coarse grid correction scheme and procedure for smoothing of 
coarse grid corrections, found crucial in the present work will be emphasized. In this 
paper we consider a Runge-Kutta scheme [6] as the smoother for the multigrid method. 
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Central differences for spatial approximations are augmented by an artificial viscosity 
based on TVD principles [7], Several changes are made to the numerical algorithm so 
that a converged solution can be obtained for high-speed flows. We initially describe 
e Runge-Kutta method for the central-difference scheme with numerical viscosity. We 
finally present some examples to demonstrate our conclusions. 


Basic Scheme 

The basic elements of the scalar dissipation model considered in this paper were 
rst intro uced by Jameson, Schmidt, and Turkel [6] in conjunction with Runge-Kutta 
explicit schemes. The spatial discretization is based on central differences with an addi- 
jonal artificial viscosity. This algorithm has been used by many investigators to solve 
e Euler equations numerically for a wide range of fluid dynamic applications. The 

/ a nn tyP u 6 ° f mI discretization has be en applied to alternating direction implicit 
( ) schemes [8] and LU factored implicit schemes [9], In this section the basic scheme 

is briefly reviewed. 

Consider the Euler equations in the form 


rrt i- Jx "h g y 


where the four-component vector of conserved variables 

W = { p pu pv pE ] T , 


( 2 ) 


and /, g are the corresponding flux vectors. The quantity p is the density, u and u 
are the Cartesian velocity components, and E is the specific total internal energy The 
independent variables are time < and Cartesian coordinates (*,y). If (1) is transformed 
to arbitrary curvilinear coordinates = {( x ,y) and r, = t}(x,y), then we obtain 

(J-'W) t + F € + G v = 0, 

where is the inverse transformation Jacobian, and 

^ = fVv ~ 9 x v> G = gx ^ — fy 

In a cell-centered, finite- volume method, (1) is integrated over an elemental volume in 
the discretized computational domain, and J” 1 is identified as the volume of the cell, 
liquation (3) can also be written as 

J~ x W t + AWi + BW V = 0, 

where A and B are the flux Jacobian matrices defined by A = dF/dW and B = dGjdW 

To advance the scheme in time we use a multistage scheme. A typical step of a 
Runge-Kutta approximation to (3) is P 

W(*> = w<°> - a k j± [D^-V + D V G w - AD) , (4) 

where D ( and A, are spatial differencing operators, and AD represents the artificial 
issipation terms. The derivatives of the fluxes are approximated by central differences. 
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The dissipation terms are a blending of second and fourth differences. That 



AD = ( 

D\ + D\ - D\ - d;) W, 

(5) 

where 

D'\W = \ 


(6) 


D\W = V* 

A<V<A( ] Wu ' 

(7) 


and At , V< are the standard forward and backward difference operators, respectively, 
associated with the £ direction. The variable scaling factor A is chosen as 


hi+y = \ ^ 

where A c is proportional to the spectral radius of the matrix A. The coefficients e( 2 > and 
e {*) are adapted to the flow and are defined as follows: 

(9) 


\Pi+l,j ~ 2 Pi.j + Pjzh 


Vij - 

\pi+ ij + 2 pij + Pt-ij | 

M) _ 

max [ 0 , (k (4) - 

and the quantities and < 


( 10 ) 

( 11 ) 


r 1 ' * • *i 

The operators for the rj direction are defined in a similar manner. 

In this paper we will also consider a matrix form of the dissipation model just de- 
scribed. The model of (5)-(ll) is characterized as a scalar formulation, since the dis- 
sipation for each discrete conservation equation is scaled by the same eigenvalue. As 
discussed in [7], a matrix form is obtained by replacing A with a matrix, so that eac 
equation will be scaled by its corresponding eigenvalue. That is, m the £ direction, the 
\A\ is substituted for the eigenvalue scaling factor, A, in (6) and (7). For the r) direction, 
£ and \A\ are replaced by t) and |J?|, respectively. A convenient form for the matrix \A\ 

is defined in the following way. Let 


A = Diag [Aj A 2 A 3 A 3 ] 


with 


a\ c, 


A 2 — q ~ V 0 ! + 0-2 C i A 3 — Qi 


Ai = q+ \[a 5 + 

a 1 = J~ 1 £ x , a 2 = q = a x u + a 2 v. 

and c representing the speed of sound. Then, 


\A\ = | A 3 \I + f 


|Ai | + | A 2 I 


I'M 


7 - 1 


E 1 + 


1 


a\ + a 2 


+ 


|Ai| — | A 2 1 


, \j a \ + a 2 


[Ej + (7 — 1)^4] 5 


( 12 ) 
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where 


E 1 = 


4> 

U(f> 

V(f> 

Hcj> 


—u 

-u 2 

—uv 

-uH 


1 

u 

— V 2 V 

-vH H 


—v 

—uv 

,2 


E 2 = 


0 


E,= 


0 

-a x q a\ 
o. 2 q a 2 a x 
-q 2 qa x 

a i 

Udi 

va i 


0 0 
(1x0,2 0 

0 
0 


aj 

qa 2 


-q fli a 2 o 

—uq ua x ua 2 0 
—vq va x va 2 0 
L —Hq Hai Ha 2 0 


E, = 


0 

ax<f> 


0 

—a x u 


0 

—a x v 
—a 2 v 
—qv 


0 
a l 
«2 
q J 


a 2 <t> —a 2 u 
q<f> —qu 

H is the total enthalpy, and <f> = ( u 2 + v 2 )/2. Notice the special form of \A\, where each 
row of Ej is either a scalar times the first row or a scalar times the second row when the 
rst row contains only zeros. Due to this special form for any A,, A 2 , andA 3 , an arbitrary 
vector * can be multiplied by \A\ very quickly. That is, one calculates | A j+ ±\ (u J+1 - Uj ) 

directly, rather than calculate \A J+ i_\ and multiply a matrix times a vector' The matrix 
\B\ is computed in the same way as \A\ by simply replacing £ with r). 

In practice one cannot choose A a ,A 2 ,A 3 as given above. Near stagnation points A 3 
approaches zero while near sonic lines Ai or A 2 approach zero. A zero artificial viscosity 
would create numerical difficulties. To prevent such problems, these values are limited 

clS 

|Ai| — maxdAil, V n p(A)), p(A) — |^| -f- c\J~a\ + a|, (13) 

|A 2 | = max(|A 2 |, V n p{A)), |A 3 | = max(|A 3 |, V e p(A)), ( 14 ) 

where the linear eigenvalue A 3 can be limited differently than the nonlinear eigenvalues 
The parameters 14 and V t are determined numerically, and the value used here is 0 25 
The second-difference term in these dissipation models is nonlinear. Its purpose is to 
mtroduce an entropy-like condition and to suppress oscillations in the neighborhood of 
s ocks This term is small in the smooth portion of the flow field. The fourth-difference 
dissipation term is basically linear and is included to damp high-frequency modes and 
a low the scheme to approach a steady state. Only this term affects the linear stability 
o e scheme. Near shocks it is reduced to zero. For high speed flows the switch (10) is 

no very good and does not allow the multigrid to converge. Instead we consider a TVD 
variation of the switch [7] given by 


\Pi±hL ~ + Pi-i.jl 


l/i t j = 1 — 1 _ — 

PiJ I I Pi,j ~ Pi—\,J I “f € 


.(2) 


= 1 / 2 . 


(15) 


With this change and the factor 1/2 in front of the second-difference dissipation term 
the scalar equation becomes first-order upwind near shocks. In the case of the original 
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v we find that v ~ .05 near shock waves in transonic flows. The parameter e must be 
chosen carefully to prevent the switch from being activated by noise. In fact we found 
it useful to take an average of the two versions for v. Hence, we use 

lP»+i,j ~ ~b 1 ,j I ^ (16) 

_ (1 - e) * (|p,+i,j - Pi,)\ + | Pi,j ~ P«-i,jl) + c * (Pi+hj + 2 pij + Pi-i, j) 

with £= 1/2 a reasonable compromise. We now no longer have a free parameter for the 
second-difference dissipation. 

Several other changes were made to the scheme in addition to the change to a TVD 
switch. In the original algorithm the artificial viscosity for the energy equation was 
based on the total enthalpy rather than the total internal energy. For high speed flows 
we base the artificial viscosity on the total internal energy so that in each equation the 
basic dependent variable is also used in the artificial viscosity. This is more in line with 
upwind schemes. This has previously been used in central-difference schemes [10]. The 
algorithm no longer preserves a constant total enthalpy in the steady state (as the Euler 
equations do), but enthalpy damping is not useful for supersonic flows. In most cases the 
difference between the two approaches is small with each approach having its advantages. 
The original form seems to give slightly sharper shocks, while the other one appears to 
make the scheme more robust. 

The form of the dissipation model of the basic or driving scheme is usually modified 
for coarse grid problems in the multigrid process. A constant coefficient second-difference 
dissipation is not only less expensive computationally but also generally provides ade- 
quate smoothing properties. For high speed flows we find it necessary, as in [4], to append 
a nonlinear dissipation to the usual one. Here this nonlinear contribution depends on 
the modified switching function of (16). We also need to increase the constant coefficient 
from the standard value of 1/16 to a value of 1/4. 

In order for the scheme to be stable it is necessary to restrict the time step. For 
transonic flows it is sufficient to base this limitation on the inviscid terms except for 
extremely fine meshes. For higher speed flows we found it necessary to include a viscous 
correction to the time step restriction even for crude meshes. We shall thus develop a 
sufficient condition for stability for the thin- layer equations. Hence, we introduce body 
fitted coordinates and then ignore all second derivatives except for the rjTj derivative. 
One then obtains the linearized equation 

W t + AWt + BW V = CW vn , (17) 


where again W is the vector of conserved variables, and the tilde indicates that a matrix 
is multiplied by the transformation Jacobian J. The matrices of (17) are somewhat 
complicated. By transforming to nonconservative variables the matrices are greatly 
simplified, and they have the same eigenvalues. Moreover, IF, A, /?, and C are replaced 

by 

W * = [ p u v p ] r , (18) 


= [P 

A* = MAM~\ B* = 

where the matrix M is defined according to 


C* = A/CM" 


(19) 


ow* 

da ~ M da ’ 


( 20 ) 
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and a is any independent variable. Abarbanel and Gottlieb [11] have shown that one 
can simultaneously symmetrize all these matrices. They symmetrize the matrices with 
the similarity transformation determined by 


S p = 


& 0 0 
0 1 0 
0 0 1 


0 

0 

0 


0 0 ^ pc 


L yfv 


( 21 ) 


S. 


- 1 


C 

y/yp 

0 

0 


( 22 ) 


0 0 0 

1 0 0 

0 1 0 

. 0 0 -A-- 

Using this transformation, the triangle inequality, and the condition for symmetric ma- 
trices that the spectral radius is equal to the norm, one can easily show that a sufficient 
stability condition is 

111 1 

(23) 


1 1 1 


At ~ A t f 


At r , 


viscous 


The quantity 1/A is bounded by the spectral radius of S 1 A*S P given by 


^ = !«£* + V&l + C \Jtl + £y, 

and similarly the quantity 1/A is bounded by 


K = I VT)y + UTf x \ + Cy/lfc + ^2. 

The matrix S^C'Sp is given by 

£_ 

3p 

and its spectral radius is given by 


0 0 

0 + 3 Tfl r) x t) y 0 

0 VxVy 47/2-1-3^2 o 

0 0 0 U(ril + ril 


(24) 


viscous 


^(Vl + V 2 y)max[^4). 


In general the first term in the maximum will be the larger. Hence, we can replace 
A v i S cou S >n (23) by its upper bound. So the actual time step (A t act ) is determined as 


follows: 


A tad < N 


IP 


1-1 


*< + *, + + vl) 


(25) 


where N is taken to be the allowable CFL number, and the constant d is 4. In the 
case of steady flows, one can advance the solution at each grid point with the time step 
determined from this estimate. This type of time stepping provides a preconditioning of 
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the matrix for the system of difference equations. The preconditioning relaxes stiffness 

due to variations in local flow properties. _ . , , , 

For all flow calculations in this paper a five stage Runge-Kutta scheme with a weigh e 
evaluation, as detailed in [12], of the dissipation terms on the first, third, and fifth stages 
is used. As described above the time step is reduced in the boundary layer by including 
the viscous contribution to the time step. In addition, the time step is reduced near 
shocks by including a term that depends on v i<3 . The reduction is constructed so t a 
there is a CFL number of 1 when v = 1. It serves to reduce the magnitude of the change 
in the solution near the shock wave, which exhibits strong nonlinear behavior. 


Implicit Residual Smoothing 

Implicit residual smoothing of the residuals is used to extend the stability range o 
the basic time-stepping scheme. For two-dimensional problems, the residual smoothing 

can be applied in the form 

(/ - <)(/ - (26) 
where the residual is defined by 


= ttm - AD< m >] , rn = 1,5 (27) 

and computed in the Runge-Kutta stage m, and AD (m) is the total artificial dissipation 
at stage m, and Hff is the final residual at stage m after the sequence of smoothings 
in the £ and tj directions. The difference operators C c and C D are associated with the 
convection and physical diffusion terms. To derive the maximum stability extension for 
the hyperbolic problem, the implicit procedure is applied after each stage of the Runge- 
Kutta scheme. The coefficients 0 € and (3 r , are variable and functions of the spectral radii 
A $ and A,. They can be written as follows: 


h 


max 


1 

4 


a * — y _ i 

N* 1 + 0 r vi ) 



Pv 


max 


'(N 

1 ^ 

\N- 

1 + ) 



(28) 


where the ratio r vi = A„/A 4 , and the quantity N/N * is the ratio of the CFL number 
of the smoothed scheme to that of the basic explicit scheme (usually having a value of 
2). In hypersonic flow applications we found it necessary for N* to be 3.25, rather than 
the value of 3.75 used for transonic computations. From a linear stability analysis, the 
scheme with these coefficients is stable for all mesh cell aspect ratios when the parameter 
0 ~ .125 and N/N * is sufficiently large. The practical limitation on the Courant number 
is due to the requirement for effective high frequency damping. For large N/N * the 
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high frequency damping of the scheme vanishes. The variable coefficients are functions 
of the local mesh cell aspect ratio, and thus the smoothing process is not activated 
in a coordinate direction where it is not needed. This is important for best possible 
convergence. For further discussion of implicit residual smoothing see [13]. 

Multigrid Method 

As indicated earlier the salient features of the multigrid method considered here are 
fairly standard. Moreover, we apply the Full Approximation Storage (FAS) scheme of 
Brandt [14] to define the equivalent fine grid problem on a coarse grid. Coarser meshes 
are obtained by eliminating every other mesh line in each coordinate direction. The 
grid transfer operators for the solution, residual, and coarse grid corrections are those 
introduced by Jameson [15]. In particular, on the auxiliary meshes, the solution is 
initialized as 


u/(°) _ Wh 

2h n 2h ' 


(29) 


where the subscript refers to the mesh spacing value, the sum is over the four fine grid 
cells that compose the 2 h grid cell, and 0 is a cell volume. This rule conserves mass, 
momentum, and energy. On a coarse grid, a forcing function P is added to the governing 
discrete equations in order to impose the fine grid approximation. After the initialization 
of the coarse grid solution, this function is computed as follows: 


P2h = E R l(W h )-R2k(wM), (30 ) 

where R h {W h ) = C h W h . Then, the time-stepping scheme on the (m + l) st stage becomes 

W ^" = - <w. ~[R,h(W^) + (31) 

We can also define a new value R m for the residual as 


R-2h — ^2 h(W 2 h) + P 2 h, (32) 

collect this value, restrict the solution W 2h to the next coarser grid, and repeat the 
process. The corrections computed on a coarse grid are transferred back to a finer grid 
with bilinear interpolation. In order to execute the multigrid strategy we employ a fixed 
W-type cycle. To provide a well conditioned starting solution for the fine mesh a Full 
Multigrid (FMG) method is used. The FMG is analogous to grid sequencing, except 
multigrid cycles are performed on each coarse grid. 

Some of the additional elements of the multigrid method are not necessarily standard 
A smoothing of the coarse grid corrections being transferred to the finest grid was found 
to be beneficial in transonic computations [12]. The smoothing was accomplished with 
the implicit residual smoothing mentioned previously and a constant coefficient f3 « 0.1. 
This smoothing of the residuals on the way to finer meshes is crucial for the convergence 
of the multigrid for hypersonic flows. Such a process acts to reduce high frequency 
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oscillations caused by the interpolation. Hence, it becomes especially important near 
strong shocks, where nonphysical upstream influence can occur. Another important 
element for high Mach number (M > 10) flows is the coarse grid correction scheme. 
That is, the physical viscous terms should also be computed on the coarse meshes. 

Boundary Conditions and Initialization 

At a solid surface (wall) boundary the no-slip condition is enforced. The wall pressure 
is set to the value at the first interior solution point, and thus, a reduced normal mo- 
mentum equation is satisfied. The wall temperature ( 7 iu ) is specified. In a finite-volume 
formulation, this amounts to treating the Cartesian velocity components and the tem- 
perature difference T — T w as antisymmetric functions with respect to the wall. For each 
of the physical problems considered the Mach number at the inflow boundary exceeds 
1.0. Consequently, the dependent variables are specified at this boundary according to 
the flow conditions. At any outflow boundary, we apply simple extrapolation of the 
components of the solution vector. In general, for hypersonic flows numerical difficulties 
are experienced at the start of a calculation if the discrete flow field is initialized with 
free-stream conditions. To avoid these difficulties we apply the following procedure. The 
Mach number of the flow is set to a lower value (i.e., 2.0) than the required one. In 
addition, the wall temperature T w is set to the free-stream value. Then the Mach num- 
ber and T w are gradually increased over a few hundred time steps until the desired flow 
conditions are obtained. This Mach number ramping is only done on the coarsest mesh 
in the FMG sequence. 

Results 

We consider two-dimensional (2-D) hypersonic laminar flow over two different geome- 
tries in order to evaluate the present multigrid method. The first geometry is a circular 
cylinder. For this case the free-stream Mach number (M^) is 6.5, and the Reynolds 
number (Rep) based on the cylinder diameter D (0.82 meters) is 1.04 x 10 s . The free- 
stream temperature (Too) is 202° Kelvin, and the wall temperature is specified at 294° 
Kelvin. This represents a fairly cold wall condition relative to the temperature after the 
normal portion of the bow shock. Computed surface pressures and heat transfer rates 
are compared with the experimental data of Wieting [16]. The second geometry is the 
15 degree compression ramp tested by Holden and Moselle [17]. For this flow problem 
the free-stream Mach number is 14.1, and the Reynolds number based on a reference 
length L (0.44 meters) is 1.04 x 10 s . The length of the flat plate preceding the ramp is 
L. The Too is 89° Kelvin, and the wall temperature is 296° Kelvin. Surface distributions 
of pressure coefficient (c p ), skin-friction coefficent ( cj ), and heat transfer coefficient (c/,) 
calculated with the present multigrid method are compared with experimental data of 
[17]. These coefficients have the standard definitions. In all calculations for both cases 
we assume that the working fluid (air) is thermally and calorically perfect. Sutherland’s 
law is used to determine the molecular viscosity. 

Cylinder Flow 

The 2-D cylinder flow was computed on a 64 x 64 grid and a 128 x 128 grid. In 
figure 1 the 64 x 64 grid, which is a proper subset of the 128 x 128 grid, is shown. 
For both grids the circumferential spacing is uniform. In the normal direction on the 
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centerline, the mesh is clustered at the surface with minimum spacings of approximately 
4 x 10 -4 D and 2 x 10 -4 D for the two meshes. At the circumferential angle ( 6 ) of -90 
degrees or +90 degrees, the normal mesh spacings are increased by nearly 60 percent 
of the centerline values. This was done to accomodate the boundary-layer growth as 
well as the resolution of the inviscid flow region. As evident from figure 2, the normal 
spacing through the shock region is uniform. Convergence histories, which define the 
variation of the error with multigrid cycles, corresponding to both grids are displayed 
in figure 3. The error is measured as an rms value of the residual for the continuity 
equation. In figure 3 one observes three out of the four levels of refinement in the FMG 
procedure. The first level, which requires only a couple of CPU seconds, is used for 
the initialization (Mach number ramping), and thus it is not shown. There are three 
grids on both the third and the fourth levels. On the 128 x 128 mesh the residual is 
reduced nearly 6 orders of magnitude in 300 cycles. This requires about 5 minutes of 
CPU time on a Cray YMP. It should be emphasized that for engineering accuracy (he., 
residual reduced by 3 orders) the finest mesh calculation required about 2 minutes. Note 
that when engineering accuracy is achieved, there is no appreciable improvement in the 
viscous solution accuracy by further residual reduction. 

In figure 4 the computed surface distributions of pressure and heat transfer rate are 
presented. There is very good agreement between the predictions with the 128 x 128 grid 
and the experimental data. For the 64 x 64 grid the scaled heat transfer rate ( Q/Q re f ) 
is overpredicted for 0 < -25 degrees and 0 > +25 degrees. This indicates that opening 
the normal mesh spacing adjacent to the surface produced a spacing too large for the 64 
x 64 grid, when only a first-order approximation is used for the temperature derivative. 
The Mach number and pressure contours for the two calculations are shown in figures 5 
and 6, respectively. 1 he smoothness of the contours is evident, and the improved shock 
resolution with mesh refinement is readily seen. In addition, one can notice that the 
boundary layer is extremely thin for this case. 

Compression Ramp Flow 

The 2-D compression ramp flow was computed on grids consisting of 56 x 64 (number 
of streamwise cells x number of normal cells) and 112 x 128 cells. Figure 7 depicts the 
56 x 64 grid. There is streamwise clustering at the leading edge of the flat plate and 
at the start of 15 degree ramp. The minimum spacing is approximately 5.8 x 10“ 3 L. 
Again, to resolve the boundary layer, the mesh is clustered in the normal direction near 
the surface. The normal spacings for the coarse and fine grids are about 4.6 x 10" 4 L 
and 2.3 x 10 4 L, respectively. In figure 8 the convergence histories for this case with 
the scalar and matrix forms of the dissipation model are presented. The convergence 
rates with both forms are quite good on the 56 x 64 grid, allowing the residual to be 
decreased about 6 orders of magnitude in 300 cycles. The average rate of reduction of 
the residual with the scalar model (0.954) is slightly faster. There is some deterioration 
in the lates with mesh refinement. This slowdown with mesh refinement is also observed 
for transonic computations. Although the average rates of residual decay using the two 
dissipation forms is essentially the same (0.965), the asymptotic rate is faster using the 
scalai model. The residual is reduced just about 5 orders in 300 multigrid cycles with 
the scalai formulation. It should be pointed out that in the multigrid calculation with 
the scalar model four grids were applied. Only three grids were used in conjunction 
with the matrix model, due to numerical difficulties caused by the sudden switch from 
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sics 

-r77r^ 

as a consequence of the flow being compressed. nres sure 

In figures 9-1 1 comparisons are made between the computed var.at.ons of the pressu«, 

skin-friction, and heat transfer coefficients and the correspon mg experim 
"»s the shock capturing capability of the present central- 
results calculated with the code developed by J. L. Thomas, wh.ch ,s based on the 
Riemann solver of Roe and described in [18], are also included in these figures. The 
computed distributions exhibit excellent agreement with the data m nearly all cases^ 
With the scalar dissipation model, there are differences beween the solutions on . the 56 
x 64 grid and the 112 X 128 grid. The results obtained using the matrix model for these 
twoTds almost coincide. Moreover, the solution computed with the matrix mod e. on 
a 56 x 64 mesh is comparable to the one calculated with the scalar formulation on - 
X 128 mesh Figures 12 and 13 show the pressure contours on the ramp for each of 
presemTcomput^tions. On. can clearly see the effects of dissipation and mesh size on 
the leading edge shock and the interaction region of the two shocks. 

C °rmulufrid e r^thod with central differencing has been successfully applied to the 
solution of hypersonic viscous flows. An explicit five stage Runge-Kutta scheme has been 
used as a smoother in solving the time-dependent, thin-layer Navier-Stokes equation^ 
In this paper considerable emphasis has been focussed on the dissipative character, 
of the driving scheme for the multigrid process. The presence of strong shocks has 
required the introduction of a switching function for the numerical l dissipation based 
on TVD principles. In addition, as a consequence of the strong s oc ' s ’ a non r 
coefficient, which is dependent on this switching function, has been included in the 
coarse grid dissipation formulation. This nonlinear coefficient is not needed for transonic 
computations. We have also considered both scalar and matrix forms of the dissipation 

™ Numerical solutions have been obtained for hypersonic laminar flow over a 2-D cylin- 
der and a 2-D compression ramp. The agreement between predictions and experimental 
data is quite good. Engineering accuracy has been obtained rapidly in all computa ions, 
requiring about 2 CPU minutes on the Cray YMP. 


AC ThTI 1 uthors would like to express their appreciation to Dr. J. L. Thomas of NASA 
Langley Research Center for providing numerical results obtained with the upwind 

scheme of Roe. 
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Figure 2 Blowup of 64 x 64 grid for 2-D cylinder flow 



Figure 3 Convergence histories for 2-D cylinder flow 
computations (Moo = 6.5, a = 0.0 deg, Rep = 1.04 x 10 5 ) 


1! 




(b) Heat transfer rate distribution 


Figure 4 Surface distributions of pressure and heat transfer rate 
for 2-D cylinder flow (Af* = 6.5, a = 0.0 deg, Re D = 1.04 x 10 5 ) 



(a) 64 x 64 grid 



(b) 128 x 128 grid 


Figure 5 Mach number contours, with AM = 0.25, for 2-D 
cylinder flow (Moo = 6.5, a = 0.0 deg, Rep = 1.04 x 10 5 ) 
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(a) 64 x 64 grid 



(b) 128 x 128 grid 


Figure 6 Pressure contours for 2-D cylinder flow 
(Mo© = 6.5, a — 0.0 deg, Rep = 1.04 xlO 5 ) 
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Figure 7 A 56 x 64 grid for 2-D compression ramp flow 


19 




Figure 8 Convergence histories for 2— D compression 
ramp flow (Afoo = 14.1, Rei = 1.04 xlO 5 , 15 degree ramp) 
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Log(100 x Cp/2) 




(b) Matrix dissipation 


Figure 9 Pressure coefficient distributions for 2-D compression 
ramp flow (Moo = 14.1, Re L = 1-04 xlO 5 , 15 degree ramp) 


21 



1 00 x Cf 1 00 x Cf/2 



(a) Scalar dissipation 



(b) Matrix dissipation 

Figure 10 Skin-friction coefficient distributions for 2-D compression 
ramp flow (M^ = 14.1, Re L = 1.04 xlO 5 , 15 degree ramp) 



(a) Scalar dissipation 



Figure 1 1 Heat transfer coefficient distributions for 2-D compression 
ramp flow (Moo = 14.1, Rei = 1.04 xlO 5 , 15 degree ramp) 
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(b) Matrix dissipation 


Figure 12 Pressure contours for 2-D compression ramp flow on 
56 x 64 grid (M^ = 14.1, Rei = 1.04x10 s , 15 degree ramp) 




1.00 r 



(a) Scalar dissipation 



(b) Matrix dissipation 

Figure 13 Pressure contours for 2-D compression ramp flow on 
112 x 128 grid (Moo = 14.1, Re L = 1.04xl0 5 , 15 degree ramp) 
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